load packages and functions
source('mytheme.R')
# model version
modelversion = 'log_rstan'
rstanmodelPath = 'modelrlt'
modelPath = paste0(rstanmodelPath, '/models/', modelversion)
Merge the model Result data
load data
AllExpData = read.csv(paste0("../data/AllValidData.csv"))
dur <- sort(unique(AllExpData$curDur))
AllExpData$WMSize <- factor(AllExpData$WMSize, labels = c("low", "medium", "high"))
# 1: 500ms, 2: 2500, 3 : 2000ms the mean reaction time of WM test
AllExpData$gap <- factor(AllExpData$gap, labels = c('short','long', 'not sure'))
AllExpData[which(AllExpData$Exp == 'Exp4'),"Exp"] = "Exp4a"
AllExpData[which(AllExpData$Exp == 'Exp5'),"Exp"] = "Exp4b"
Corrct rate
WMCrr2 <- ggplot(meanForPlot, aes(WMSize, mean_WMCrr, ymin = mean_WMCrr - se_WMCrr, ymax = mean_WMCrr + se_WMCrr,
group =Exp, color = Exp, fill = Exp))+
geom_line(stat = "identity",position = position_dodge(width = 0.2))+
geom_point(stat = "identity",position = position_dodge(width = 0.2))+
geom_errorbar(width=.2, position = position_dodge(width = 0.2)) +
coord_cartesian(ylim = c(0.5, 1)) +
colorSet5+
labs(x = "Memory load", y = "Mean accuracy in WM task") +
theme_new
WMCrr2
ggsave(paste0(getwd(), "/figures/WMCrr2.png"), WMCrr2, width = 4, height = 4)
### generate WM correct rates
AllExpData$WMCrr <- AllExpData$TPresent == AllExpData$WMRP
m_wmp<- dplyr::group_by(AllExpData, Exp, WMSize, NSub) %>%
dplyr::summarize(m_WMCrr = mean(WMCrr), n =n(), se_WMCrr = sd(WMCrr)/sqrt(n-1))
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
comparison between Exp.4a and Exp.4b)
correct rate
#plot WM correct rates
dplyr::group_by(AllExpData%>%filter(Exp %in% c('Exp4a', 'Exp4b')), Exp, WMSize, NSub, gap) %>%
dplyr::summarize(m_WMCrr = mean(WMCrr), n = n(), se_WMCrr = sd(WMCrr)/sqrt(n-1)) -> correctrate_gap
## `summarise()` has grouped output by 'Exp', 'WMSize', 'NSub'. You can override
## using the `.groups` argument.
write.csv(correctrate_gap, paste0(modelPath, '/rlt/correctrate_gap.csv'))
correctrate_gap%>%
dplyr::group_by(Exp, WMSize, gap)%>%
dplyr::summarize( n = n(),
mean_WMCrr = mean(m_WMCrr), se_WMCrr = sd(m_WMCrr)/sqrt(n-1) ) -> meanForPlot_gap
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
WMCrr_gap <- ggplot(meanForPlot_gap, aes(WMSize, mean_WMCrr, ymin = mean_WMCrr - se_WMCrr, ymax = mean_WMCrr + se_WMCrr, group =interaction(Exp,gap), color = gap, shape = Exp, linetype = Exp))+
geom_line(stat = "identity",position = position_dodge(width = 0.2))+
geom_point(stat = "identity",position = position_dodge(width = 0.2))+
geom_errorbar(width=.2, position = position_dodge(width = 0.2)) +
coord_cartesian(ylim = c(0.5, 1)) +
colorSet5+
labs(x = "Memory load", y = "Mean accuracy in WM task") +theme_new
WMCrr_gap
ggsave(paste0(getwd(), "/figures/WMCrr_gap.png"), WMCrr_gap, width = 4, height = 4)
observed RP
AllExpData%>% filter(Exp %in% c('Exp4a', 'Exp4b')) %>%
dplyr::group_by(Exp, curDur, WMSize, NSub, gap) %>%
dplyr::summarize(n = n(),
m_repDur = mean(repDur),
se_repDur = sd(repDur)/sqrt(n-1)) ->RP_obs_gap_bias
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize', 'NSub'. You can
## override using the `.groups` argument.
RP_obs_gap_bias$rep_bias = RP_obs_gap_bias$m_repDur - RP_obs_gap_bias$curDur
ezANOVA(data = RP_obs_gap_bias%>%filter(Exp =='Exp4b'), dv= rep_bias, wid=NSub, within= .(curDur, gap, WMSize) )
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: "curDur" will be treated as numeric.
## Warning: You have removed one or more levels from variable "gap". Refactoring
## for ANOVA.
## Warning: There is at least one numeric within variable, therefore aov() will be
## used for computation and no assumption checks will be obtained.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 1 curDur 1 15 78.8694018 2.325297e-07 * 8.072483e-01
## 2 gap 1 15 0.3640022 5.553058e-01 5.913116e-04
## 3 WMSize 2 30 9.8759805 5.066007e-04 * 6.980625e-02
## 4 curDur:gap 1 15 10.6242360 5.279692e-03 * 1.079097e-02
## 5 curDur:WMSize 2 30 13.1667406 7.858144e-05 * 2.343741e-02
## 6 gap:WMSize 2 30 3.1032330 5.956988e-02 2.529934e-03
## 7 curDur:gap:WMSize 2 30 0.1127064 8.937894e-01 7.606329e-05
RP_obs_gap <- ggplot( data = RP_obs_gap_bias%>%
dplyr::group_by(Exp, curDur, WMSize, gap) %>%
dplyr::summarize(m_m_repDur = mean(m_repDur),
m_se_repDur = mean(se_repDur)), aes(x = curDur, y = m_m_repDur, color=as.factor(WMSize), shape = gap)) +
geom_point()+
geom_line()+
geom_abline(slope=1, intercept=0)+
facet_grid(cols = vars(Exp)) +
labs(x="Sample intervals (s)", y="Reproduction (s)", shape=" ", color = "Memory Load")+
theme_new+colorSet3+
scale_x_continuous(breaks=seq(0, 1.6, 0.4))+ theme(legend.position="top")
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize'. You can override
## using the `.groups` argument.
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_obs.png"), RP_obs_gap, width = 6, height = 4)
RP_obs_gap
central tendency effect (slope) and IPs (linear regression)
#Observed Indifference Point
obs_model <- function(df) {
lm(repDur ~ curDur, data = df)
}
obs_InP <- AllExpData %>% filter(Exp %in% c('Exp4a', 'Exp4b')) %>%
dplyr::group_by(NSub, Exp, WMSize, gap) %>% nest() %>%
mutate(model = map(data, obs_model)) %>% # linear regression
mutate(slope = map(model, broom::tidy)) %>% # get estimates
unnest(slope, .drop = TRUE) %>% # remove raw data
select(-std.error,-statistic, -p.value) %>% # remove unnessary columns
spread(term, estimate) %>% # spread stimates
dplyr::rename(Intercept = `(Intercept)`, slope = curDur) # rename columns
## Warning: The `.drop` argument of `unnest()` is deprecated as of tidyr 1.0.0.
## All list-columns are now preserved.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
obs_InP$model = NULL
obs_InP$data = NULL
obs_InP$inP = obs_InP$Intercept /(1-obs_InP$slope)
write.csv(obs_InP, paste0(modelPath, '/rlt/obs_InP_gap.csv'))
mobs_InP = obs_InP %>% group_by(Exp, WMSize, gap)%>%
dplyr::summarise(n=n(),
m_Intercept = mean(Intercept),
se_Intercept= sd(Intercept)/sqrt(n-1),
m_inP = mean(inP),
se_inP = sd(inP)/sqrt(n-1),
m_slope = mean(slope),
se_slope = sd(slope)/sqrt(n-1))
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
plt_InP_linear_gap<- ggplot(data = mobs_InP, aes(x=WMSize, y=m_inP, group = gap, color = gap))+
geom_line(stat = "identity",position = position_dodge(width = 0.2))+
geom_point(stat = "identity",position = position_dodge(width = 0.2))+
geom_errorbar(width=.3, aes(ymin = m_inP - se_inP, ymax = m_inP + se_inP), position = position_dodge(width = 0.2)) +theme_new+
labs(colour = "Gap")+colorSet3+
facet_wrap(~Exp)+
xlab(' ')+ylab("indifference point (s)")+guides(shape="none")+
theme(legend.position = "top")
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_InP_linear_gap.png"), plt_InP_linear_gap, width = 5, height = 5)
plt_InP_linear_gap
ezANOVA(data = obs_InP%>%filter(Exp =='Exp4b'), dv= inP, wid=NSub, within= .(gap, WMSize) )
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: You have removed one or more levels from variable "gap". Refactoring
## for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 gap 1 15 2.942083 0.106876928 0.006249289
## 3 WMSize 2 30 7.307023 0.002598632 * 0.069252661
## 4 gap:WMSize 2 30 1.289560 0.290220752 0.006533781
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 3 WMSize 0.3727819 0.001000422 *
## 4 gap:WMSize 0.5400426 0.013396643 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 3 WMSize 0.6145458 0.01069076 * 0.6414052 0.009678944 *
## 4 gap:WMSize 0.6849515 0.28394478 0.7306789 0.285474638
plt_RP_slope_linear_gap<- ggplot(data = mobs_InP, aes(x= WMSize, y=m_slope, group = gap,color = gap, shape = gap))+
geom_line(stat = "identity",position = position_dodge(width = 0.2))+
geom_point(stat = "identity",position = position_dodge(width = 0.2))+
geom_errorbar(width=.3, aes(ymin = m_slope - se_slope, ymax = m_slope + se_slope), position = position_dodge(width = 0.2)) +theme_new+
facet_wrap(~Exp)+
labs(colour = "Gap")+colorSet3+
xlab(' ')+ylab("Slope of reproduction (s)")+
theme(legend.position = "top")
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_RP_slope_linear_gap.png"), plt_RP_slope_linear_gap, width = 5, height = 5)
plt_RP_slope_linear_gap
ezANOVA(data = obs_InP%>%filter(Exp =='Exp4b'), dv= slope, wid=NSub, within= .(gap, WMSize) )
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: You have removed one or more levels from variable "gap". Refactoring
## for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 gap 1 15 10.6363088 5.259698e-03 * 1.264088e-02
## 3 WMSize 2 30 13.1535570 7.913522e-05 * 2.751879e-02
## 4 gap:WMSize 2 30 0.1142223 8.924457e-01 9.048962e-05
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 3 WMSize 0.8539487 0.3311481
## 4 gap:WMSize 0.9665181 0.7878977
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 3 WMSize 0.8725613 0.000189114 * 0.9778272 9.205975e-05 *
## 4 gap:WMSize 0.9676028 0.886521995 1.1084480 8.924457e-01
load model results
check model parameters
### Average Parameters
mm_Baypar <- dplyr::group_by(Bayparlist, Exp) %>%
dplyr::summarize( m_sig_s2 = mean(sig_s2),
m_sig_pr2_log = mean(sig_pr2_log),
m_ks= mean(ks),
m_kr = mean(kr),
m_ls = mean(ls),
m_mu_pr_log= mean(mu_pr_log),
n= n(),
se_sig_s2 = sd(sig_s2)/sqrt(n-1),
se_sig_mn2 = sd(sig_mn2)/sqrt(n-1),
se_sig_pr2_log = sd(sig_pr2_log)/sqrt(n-1),
se_ks = sd(ks)/sqrt(n-1),
se_kr = sd(kr)/sqrt(n-1),
se_ls =sd(ls)/sqrt(n-1),
se_mu_pr_log = sd(mu_pr_log)/sqrt(n-1))
mm_Baypar
Average Parameters
mBaypar_sub <- dplyr::group_by(Bayparlist, Exp, NSub) %>%
dplyr::summarize( m_sig_s2 = mean(sig_s2),
m_sig_pr2_log = mean(sig_pr2_log),
m_ks= mean(ks),
m_kr = mean(kr),
m_ls = mean(ls),
m_mu_pr_log= mean(mu_pr_log),
n= n(),
se_sig_s2 = sd(sig_s2)/sqrt(n-1),
se_sig_mn2 = sd(sig_mn2)/sqrt(n-1),
se_sig_pr2_log = sd(sig_pr2_log)/sqrt(n-1),
se_ks = sd(ks)/sqrt(n-1),
se_kr = sd(kr)/sqrt(n-1),
se_ls =sd(ls)/sqrt(n-1),
se_mu_pr_log = sd(mu_pr_log)/sqrt(n-1))
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
fig_ks_subj = ggplot(mBaypar_sub, aes(Exp, m_ks, color = Exp)) +
geom_boxplot(position = position_dodge()) +
geom_jitter(shape=16, position=position_jitter(0.2))+
#facet_wrap(~Exp)+
theme_new+ colorSet5 +
theme(strip.background = element_blank()) + # remove subtitle background
labs(x = ' ', y = 'scale factor Ks')+
theme(legend.position='bottom')
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig_ks_subj.png"), fig_ks_subj, width = 6, height = 6)
fig_ks_subj
fig_kr_subj = ggplot(mBaypar_sub, aes(Exp, m_kr, color = Exp)) +
geom_boxplot(position = position_dodge()) +
geom_jitter(shape=16, position=position_jitter(0.2))+
#facet_wrap(~Exp)+
theme_new+ colorSet5 +
theme(strip.background = element_blank()) + # remove subtitle background
labs(x = ' ', y = 'scale factor Kr')+
theme(legend.position='bottom')
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig_kr_subj.png"), fig_kr_subj, width = 6, height = 6)
fig_kr_subj
fig_ls_subj = ggplot(mBaypar_sub, aes(Exp, m_ls, color = Exp)) +
geom_boxplot(position = position_dodge()) +
geom_jitter(shape=16, position=position_jitter(0.2))+
#facet_wrap(~Exp)+
theme_new+ colorSet5 +
theme(strip.background = element_blank()) + # remove subtitle background
labs(x = ' ', y = 'scale factor ls')+
theme(legend.position='bottom')
ggsave(paste0(getwd(), "/", modelPath, "/figures/mfig_ls_subj.png"), fig_ls_subj, width = 6, height = 6)
fig_ls_subj
Prediction results
#load prediction
AllDat_predY <- read.csv(paste0(modelPath, "/rlt/AllDat_predY_",modelversion,".csv"))
AllDat_predY$WMSize <- as.factor(AllDat_predY$WMSize)
levels(AllDat_predY$WMSize) = c("low", "medium", "high")
AllDat_predY$pred_Bias = AllDat_predY$mu_r - AllDat_predY$curDur
AllDat_predY$predErr = AllDat_predY$mu_r - AllDat_predY$repDur
AllDat_predY$relatErr = AllDat_predY$predErr / AllDat_predY$repDur
AllDat_predY[which(AllDat_predY$Exp == "Exp4"),"Exp"] = "Exp4a"
AllDat_predY[which(AllDat_predY$Exp == "Exp5"),"Exp"] = "Exp4b"
#calculate the mean reproduction biases for the five given intervals for all subjects
mpredY_sub <- dplyr::group_by(AllDat_predY, curDur, Exp, NSub, WMSize) %>%
dplyr::summarize(n = n(),
m_repDur = mean(repDur),
sd_repDur = sd(repDur),
m_mu_r = mean(mu_r),
m_sig_r = mean(sig_r),
m_wp = mean(wp),
se_wp = sd(wp)/sqrt(n-1),
log_lik =mean(log_lik),
cv =sd_repDur/ m_repDur,
pred_cv = mean(sig_r/mu_r),
predRP_err = mean(m_mu_r-m_repDur),
predVar_err = mean(m_sig_r-sd_repDur),
predRP_rerr = mean(abs(m_mu_r-m_repDur)/m_repDur),
predVar_rerr = mean(abs(m_sig_r-sd_repDur)/sd_repDur),
predcv_err = pred_cv-cv,
predcv_rerr = mean(abs(pred_cv-cv)/cv))
## `summarise()` has grouped output by 'curDur', 'Exp', 'NSub'. You can override
## using the `.groups` argument.
mpredY_sub_split <-split(mpredY_sub %>%select(c('Exp', 'NSub', 'WMSize', 'm_repDur', 'curDur')), mpredY_sub$curDur)
mpredY_sub_cv_split <-split(mpredY_sub %>%select(c('Exp', 'NSub', 'WMSize', 'cv', 'curDur')), mpredY_sub$curDur)
mpredY_sub_split1 <-split(mpredY_sub %>%select(c('Exp', 'NSub', 'WMSize', 'curDur' ,'m_repDur')), mpredY_sub$WMSize)
mpredY_sub_cv_split1 <-split(mpredY_sub %>%select(c('Exp', 'NSub', 'WMSize','curDur', 'cv')), mpredY_sub$WMSize)
mpredY_sub_jasp_RP = NULL
mpredY_sub_jasp_cv = NULL
mpredY_sub_jasp_RP1 = NULL
mpredY_sub_jasp_cv1 = NULL
for(i in 1: length(mpredY_sub_split)){
temp = mpredY_sub_split[[i]]
curDur = unique(temp$curDur)
temp$curDur = NULL
colnames(temp) = c('Exp', 'NSub', 'WMSize', paste0('m_repDur_', curDur))
temp_cv = mpredY_sub_cv_split[[i]]
curDur = unique(temp_cv$curDur)
temp_cv$curDur = NULL
colnames(temp_cv) = c('Exp', 'NSub', 'WMSize', paste0('cv_', curDur))
if(i == 1){
mpredY_sub_jasp_RP = temp
mpredY_sub_jasp_cv = temp_cv
}
else{
mpredY_sub_jasp_RP = left_join(mpredY_sub_jasp_RP, temp, by=c("Exp", "NSub", "WMSize"))
mpredY_sub_jasp_cv = left_join(mpredY_sub_jasp_cv, temp_cv, by=c("Exp", "NSub", "WMSize"))
}
}
for (i in 1: length(mpredY_sub_split1)){
temp1 = mpredY_sub_split1[[i]]
WMSize = unique(temp1$WMSize)
temp1$WMSize = NULL
colnames(temp1) = c('Exp', 'NSub', 'curDur', paste0('m_repDur_', WMSize))
temp_cv1 = mpredY_sub_cv_split1[[i]]
WMSize = unique(temp_cv1$WMSize)
temp_cv1$WMSize = NULL
colnames(temp_cv1) = c('Exp', 'NSub', 'curDur', paste0('cv_', WMSize))
if(i == 1){
mpredY_sub_jasp_RP1 = temp1
mpredY_sub_jasp_cv1 = temp_cv1
}
else{
mpredY_sub_jasp_RP1 = left_join(mpredY_sub_jasp_RP1, temp1, by=c("Exp", "NSub", "curDur"))
mpredY_sub_jasp_cv1 = left_join(mpredY_sub_jasp_cv1, temp_cv1, by=c("Exp", "NSub", "curDur"))
}
}
write_csv(mpredY_sub_jasp_RP, paste0(modelPath, '/rlt/RP_Bias_jasp.csv'))
write_csv(mpredY_sub_jasp_cv, paste0(modelPath, '/rlt/cv_jasp.csv'))
write_csv(mpredY_sub_jasp_RP1, paste0(modelPath, '/rlt/RP_Bias_jasp1.csv'))
write_csv(mpredY_sub_jasp_cv1, paste0(modelPath, '/rlt/cv_jasp1.csv'))
write_csv(dplyr::group_by(mpredY_sub, curDur, NSub) %>%
dplyr::summarize(m_cv = mean(cv))%>%spread(curDur, m_cv), paste0(modelPath, '/rlt/m_cv.csv'))
## `summarise()` has grouped output by 'curDur'. You can override using the
## `.groups` argument.
mpredY_sub$RP_bias = mpredY_sub$m_repDur -mpredY_sub$curDur
mpredY_sub_new <- dplyr::group_by(mpredY_sub, curDur, Exp, NSub) %>%
dplyr::summarize(m_RP_bias = mean(RP_bias))%>% spread(curDur, m_RP_bias)
## `summarise()` has grouped output by 'curDur', 'Exp'. You can override using the
## `.groups` argument.
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp1'), paste0(modelPath, '/rlt/RP_Bias_exp1.csv'))
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp2'), paste0(modelPath, '/rlt/RP_Bias_exp2.csv'))
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp3'), paste0(modelPath, '/rlt/RP_Bias_exp3.csv'))
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp4a'), paste0(modelPath, '/rlt/RP_Bias_exp4a.csv'))
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp4b'), paste0(modelPath, '/rlt/RP_Bias_exp4b.csv'))
mpredY_sub_WMsize <- dplyr::group_by(mpredY_sub, WMSize, Exp, NSub) %>%
dplyr::summarize(m_RP_bias = mean(RP_bias))%>% spread(WMSize, m_RP_bias)
## `summarise()` has grouped output by 'WMSize', 'Exp'. You can override using the
## `.groups` argument.
write_csv(mpredY_sub_WMsize%>%filter(Exp == 'Exp3'), paste0(modelPath, '/rlt/RP_Bias_WMsize_exp3.csv'))
write_csv(mpredY_sub_WMsize%>%filter(Exp == 'Exp4a'), paste0(modelPath, '/rlt/RP_Bias_WMsize_exp4a.csv'))
write_csv(mpredY_sub_WMsize%>%filter(Exp == 'Exp4b'), paste0(modelPath, '/rlt/RP_Bias_WMsize_exp4b.csv'))
#### predicted data
m_predY <- mpredY_sub%>%
dplyr::group_by(Exp, curDur, WMSize) %>%
dplyr::summarize(m_m_repDur = mean(m_repDur),
m_sd_repDur = mean(sd_repDur),
m_m_sig_r =mean(m_sig_r),
m_m_mu_r = mean(m_mu_r),
m_m_wp = mean(m_wp),
n = n(),
m_se_wp = sd(se_wp)/sqrt(n-1),
log_lik =mean(log_lik),
mpredRP_err = mean(predRP_err),
mpredVar_err = mean(predVar_err),
mpredRP_rerr = mean(predRP_rerr),
mpredVar_rerr = mean(predVar_rerr),
cv= mean(cv),
pred_cv = mean(pred_cv),
mpredcv_err = mean(predcv_err),
mpredcv_rerr = mean(predcv_rerr))
m_predY_acc = mpredY_sub%>%
dplyr::group_by(Exp) %>%
dplyr::summarize(mpred_rerr = mean(predRP_rerr)*100,
mpredVar_rerr = mean(predVar_rerr)*100,
mpredcv_rerr = mean(predcv_rerr)*100)
m_predY_acc
WAIC and LOO-CV
#extract waic and loo-cv from parameter list
m_WAIC <- dplyr::group_by(Bayparlist, Exp) %>%
dplyr::summarize(n =n(),
m_looic = mean(looic),
m_waic = mean(waic),
se_waic = sd(waic)/sqrt(n-1),
se_looic = sd(looic)/sqrt(n-1),
m_p_loo = mean(p_loo),
m_elpd_loo = mean(elpd_loo),
m_se_looic = mean(se_looic),
m_se_p_loo = mean(se_p_loo),
m_p_waic = mean(p_waic),
m_se_waic = mean(se_waic))
m_WAIC
#load test results
AllDat_newY <- read.csv(paste0(modelPath, "/rlt/AllDat_newY_",modelversion,".csv"))
AllDat_newY$WMSize <- as.factor(AllDat_newY$WMSize)
levels(AllDat_newY$WMSize) = c("low", "medium", "high")
AllDat_newY[which(AllDat_newY$Exp == "Exp4"), "Exp"] = "Exp4a"
AllDat_newY[which(AllDat_newY$Exp == "Exp5"), "Exp"] = "Exp4b"
Plot Results
RP
RP_bias_obs <- ggplot(data = AllExpData%>% filter(Exp != 'Exp5') %>%
dplyr::group_by(Exp, curDur, WMSize, NSub) %>%
dplyr::summarize(n = n(),
m_repDur = mean(repDur),
se_repDur = sd(repDur)/sqrt(n-1)) %>%
dplyr::group_by(Exp, curDur, WMSize) %>%
dplyr::summarize(m_m_repDur = mean(m_repDur),
m_se_repDur = mean(se_repDur)), aes(x = curDur, y = m_m_repDur-curDur, color=as.factor(WMSize))) +
geom_point()+
geom_line()+
#geom_errorbar(width=.2, aes(ymin = m_m_repDur-curDur - m_se_repDur, ymax = m_m_repDur -curDur + m_se_repDur)) +
geom_hline(yintercept = 0, linetype='dashed')+
facet_grid(cols = vars(Exp)) +
labs(x="Sample intervals (s)", y="Reproduction bias(s)", shape=" ", color = "Memory Load")+
theme_new+colorSet3+guides(shape="none")+
scale_x_continuous(breaks=seq(0, 1.6, 0.4))+ theme(legend.position="top")
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'Exp', 'curDur'. You can override using the `.groups` argument.
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_bias_obs.png"), RP_bias_obs, width = 6, height = 4)
RP_bias_obs
RP_bias_obs <- ggplot(data = AllExpData%>% filter(Exp != 'Exp5') %>%
dplyr::group_by(Exp, curDur, WMSize, NSub) %>%
dplyr::summarize(n = n(),
m_repDur = mean(repDur),
se_repDur = sd(repDur)/sqrt(n-1)) %>%
dplyr::group_by(Exp, curDur, WMSize) %>%
dplyr::summarize(m_m_repDur = mean(m_repDur),
m_se_repDur = mean(se_repDur)), aes(x = curDur, y = m_m_repDur-curDur, color=as.factor(WMSize))) +
geom_point()+
geom_line()+
geom_errorbar(width=.2, aes(ymin = m_m_repDur-curDur - m_se_repDur, ymax = m_m_repDur -curDur + m_se_repDur)) +
geom_hline(yintercept = 0, linetype='dashed')+
facet_grid(cols = vars(Exp)) +
labs(x="Sample intervals (s)", y="Reproduction bias(s)", shape=" ", color = "Memory Load")+
theme_new+colorSet3+guides(shape="none")+
scale_x_continuous(breaks=seq(0, 1.6, 0.4))+ theme(legend.position="top")
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'Exp', 'curDur'. You can override using the `.groups` argument.
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_bias_obs.png"), RP_bias_obs, width = 6, height = 4)
RP_bias_obs
RP <- ggplot(data = m_predY, aes(x = curDur, y = m_m_repDur, color=WMSize, shape = as.factor('Observation'))) +
geom_point(size=2, alpha = 0.5)+
geom_line(data= m_newY, aes(x=curDur, y=m_mu_r, color=WMSize)) +
geom_abline(slope=1, intercept=0)+
facet_grid(cols = vars(Exp)) +
labs(x="Sample intervals (s)", y="Reproduction (s)", shape=" ", color = "Memory Load")+
theme_new+colorSet3+guides(shape="none")+theme_new +theme(legend.position="top")
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP.png"), RP, width = 6, height = 6)
RP
RP Slope
# calculate the slope of the cv curve
RPslope_model <- function(df) {
lm(m_repDur ~ log(curDur), data = df)
}
RPslopes <- mpredY_sub %>%
dplyr::group_by(NSub, Exp, WMSize) %>% nest() %>%
mutate(model = map(data, RPslope_model)) %>%
mutate(slope = map(model, broom::tidy)) %>% # get estimates
unnest(slope, .drop = TRUE) %>% # remove raw data
select(-std.error,-statistic, -p.value) %>% # remove unnessary columns
spread(term, estimate) %>% # spread stimates
dplyr::rename(Intercept = `(Intercept)`, slope = `log(curDur)`)
RPslopes$data <- NULL
RPslopes$model <- NULL
mRPslopes <- RPslopes%>% dplyr::group_by(WMSize, Exp) %>%
dplyr::summarize(m_Intercept = mean(Intercept),
m_slope = mean(slope),
n = n(),
se_slope = sd(slope)/sqrt(n-1),
se_Intercept = sd(Intercept)/sqrt(n-1))
## `summarise()` has grouped output by 'WMSize'. You can override using the
## `.groups` argument.
plt_CVslope <- ggplot(mRPslopes, aes(Exp, m_slope, ymin = m_slope - se_slope, ymax = m_slope + se_slope, group =interaction(Exp, WMSize), color = factor(WMSize)), shape = factor(WMSize))+
geom_line(stat = "identity",position = position_dodge(width = 0.2))+
geom_point(stat = "identity",position = position_dodge(width = 0.2))+
geom_errorbar(width=.2, position = position_dodge(width = 0.2)) +
colorSet5+
labs(x = "", y = TeX("Slope of RP"), color = 'Memory Load') +
theme_new +theme(legend.position="top")
plt_CVslope
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
Anova analysis on slopes of RP
RPslopes$WMSize = as.factor(RPslopes$WMSize)
ezANOVA(data = RPslopes, dv= slope, wid=NSub, within=.(WMSize), between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 Exp 4 75 1.916369 1.164556e-01 0.08919332
## 3 WMSize 2 150 26.957878 9.954868e-11 * 0.01482364
## 4 Exp:WMSize 8 150 6.842162 1.194810e-07 * 0.01504611
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 3 WMSize 0.937441 0.09160648
## 4 Exp:WMSize 0.937441 0.09160648
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 3 WMSize 0.9411242 3.125487e-10 * 0.9645461 1.982414e-10 *
## 4 Exp:WMSize 0.9411242 2.602981e-07 * 0.9645461 1.909273e-07 *
ezANOVA(data = RPslopes %>% filter(Exp =='Exp1'), dv= slope, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 WMSize 2 30 1.628257 0.2131416 0.002398051
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 WMSize 0.8320858 0.2761702
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 WMSize 0.8562273 0.2172718 0.9557549 0.214484
ezANOVA(data = RPslopes %>% filter(Exp =='Exp2'), dv= slope, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 WMSize 2 30 21.20007 1.822755e-06 * 0.07646424
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 WMSize 0.8461525 0.3105563
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 WMSize 0.8666656 7.302883e-06 * 0.9698478 2.493639e-06 *
ezANOVA(data = RPslopes %>% filter(Exp =='Exp3'), dv= slope, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 WMSize 2 30 0.4137346 0.6648914 0.00152841
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 WMSize 0.9800567 0.8684769
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 WMSize 0.9804466 0.6610089 1.126392 0.6648914
ezANOVA(data = RPslopes %>% filter(Exp =='Exp4a'), dv= slope, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 WMSize 2 30 7.246942 0.002705914 * 0.01861911
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 WMSize 0.9730381 0.8258648
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 WMSize 0.973746 0.002974256 * 1.117022 0.002705914 *
ezANOVA(data = RPslopes %>% filter(Exp =='Exp4b'), dv= slope, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 WMSize 2 30 11.66985 0.0001782626 * 0.02763263
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 WMSize 0.8407011 0.2968187
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 WMSize 0.8625903 0.0004111601 * 0.9643404 0.0002213399 *
Anova analysis on Intercept of RP
ezANOVA(data = RPslopes, dv= Intercept, wid=NSub, within=.(WMSize), between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 Exp 4 75 1.100218 3.628149e-01 0.050380425
## 3 WMSize 2 150 4.729965 1.018536e-02 * 0.006009338
## 4 Exp:WMSize 8 150 6.814933 1.283386e-07 * 0.033669273
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 3 WMSize 0.7988765 0.0002464606 *
## 4 Exp:WMSize 0.7988765 0.0002464606 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 3 WMSize 0.8325539 1.506819e-02 * 0.8491741 1.449261e-02 *
## 4 Exp:WMSize 0.8325539 1.166472e-06 * 0.8491741 9.364511e-07 *
RP Bias
RP_bias <- ggplot(data = m_predY, aes(x = curDur, y = m_m_repDur - curDur,color=WMSize, shape = as.factor('Observation'))) +
geom_point(size=2, alpha = 0.5)+
geom_line(data= m_newY, aes(x=curDur, y=m_mu_r-curDur, color=WMSize)) +
geom_hline(yintercept = 0, linetype='dashed')+
facet_grid(cols = vars(Exp)) +
labs(x=" ", y="Reproduction bias (s)", shape=" ", color = "Memory Load")+theme_new+
colorSet3+guides(shape="none")
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_bias.png"), RP_bias, width = 6, height = 6)
RP_bias
CV
curDurItem <- unique(m_predY$curDur)
RP_CV <- ggplot(data= m_predY, aes(x=curDur, y= cv, color=WMSize, shape = as.factor('Observation'))) +
geom_point(size=2, alpha = 0.5)+
geom_line(data = m_newY, aes(x=curDur, y= m_sig_r/m_mu_r, color=WMSize)) +
facet_grid(~Exp) +
labs(x="Duration (s)", y="CV", shape=" ", color = "Memory Load")+ theme_new+
colorSet3+guides(shape="none")+theme(strip.text.x = element_blank())
RP_CV
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_CV.png"), RP_CV, width = 6, height = 6)
CV slope
# calculate the slope of the cv curve
cvSlope_model <- function(df) {
lm(log(cv_obs) ~ log(curDur), data = df)
}
mpredY <- dplyr::group_by(AllDat_predY, curDur, WMSize, Exp, NSub) %>%
dplyr::summarize(m_repDur = mean(repDur),
sd_repDur = sd(repDur),
n = n(),
sd_repDur = sd(repDur),
m_mu_r = mean(mu_r),
m_sig_r = mean(sig_r),
wp = mean(wp),
log_lik =mean(log_lik))
## `summarise()` has grouped output by 'curDur', 'WMSize', 'Exp'. You can override
## using the `.groups` argument.
mpredY$cv_obs <- mpredY$sd_repDur/mpredY$m_repDur
CVslopes <- mpredY %>%
dplyr::group_by(NSub, Exp, WMSize) %>% nest() %>% # nested data
mutate(model = map(data, cvSlope_model)) %>% # linear regression
mutate(slope = map(model, broom::tidy)) %>% # get estimates
unnest(slope, .drop = TRUE) %>% # remove raw data
select(-std.error,-statistic, -p.value) %>% # remove unnessary columns
spread(term, estimate) %>% # spread stimates
dplyr::rename(obs_cv_Intercept = `(Intercept)`, obs_cv_slope = `log(curDur)`) # rename columns
CVslopes$data = NULL
CVslopes$model = NULL
#change the table struction of slopes for spss
CVslopes_list <-split(CVslopes, CVslopes$WMSize)
CVslopes_spss = NULL
for (i in 1: length(CVslopes_list)){
temp = CVslopes_list[[i]]
WMSize = unique(temp$WMSize)
temp$WMSize = NULL
colnames(temp) = c('Exp', 'NSub', paste0('obs_cv_Intercept_',WMSize), paste0('obs_cv_slope_',WMSize))
if(i == 1)
CVslopes_spss = temp
else
CVslopes_spss = left_join(CVslopes_spss, temp, by=c("Exp", "NSub"))
}
mCVslopes <- CVslopes%>% dplyr::group_by(WMSize, Exp) %>%
dplyr::summarize(n =n(),
m_obs_cv_Intercept = mean(obs_cv_Intercept),
m_obs_cv_slope = mean(obs_cv_slope),
se_obs_CV_slope = sd(obs_cv_slope)/sqrt(n-1),
se_obs_CV_Intercept = sd(obs_cv_Intercept)/sqrt(n-1))
## `summarise()` has grouped output by 'WMSize'. You can override using the
## `.groups` argument.
plt_CVslope <- ggplot(mCVslopes, aes(Exp, m_obs_cv_slope, ymin = m_obs_cv_slope - se_obs_CV_slope, ymax = m_obs_cv_slope + se_obs_CV_slope, group =interaction(Exp, WMSize), color = factor(WMSize)), shape = factor(WMSize))+
geom_line(stat = "identity",position = position_dodge(width = 0.2))+
geom_point(stat = "identity",position = position_dodge(width = 0.2))+
geom_errorbar(width=.2, position = position_dodge(width = 0.2)) +
colorSet5+
labs(x = "", y = TeX("Slope of CV"), color = 'Memory Load') +
theme_new +theme(legend.position="top")
plt_CVslope
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_CVslope.png"), plt_CVslope, width = 4, height = 4)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
plt_CVIntecept<- ggplot(mCVslopes, aes(Exp, m_obs_cv_Intercept, ymin = m_obs_cv_Intercept - se_obs_CV_Intercept, ymax = m_obs_cv_Intercept + se_obs_CV_Intercept, group =interaction(Exp, WMSize), color = factor(WMSize)), shape = factor(WMSize))+
geom_line(stat = "identity",position = position_dodge(width = 0.2))+
geom_point(stat = "identity",position = position_dodge(width = 0.2))+
geom_errorbar(width=.2, position = position_dodge(width = 0.2)) +
colorSet5+
labs(x = "", y = TeX("Intercept of CV"), color = 'Memory Load') +
theme_new +theme(legend.position="top")
plt_CVIntecept
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_Intercept.png"), plt_CVIntecept, width = 4, height = 4)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## Figures in the MS
fig_CV_slope_Intecept<-ggarrange(plt_CVslope, plt_CVIntecept, common.legend = TRUE, ncol=2, nrow=1, labels = c("a", "b"))
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig_CV_slope_Intecept.png"), fig_CV_slope_Intecept, width = 8, height = 4)
fig_CV_slope_Intecept
## Figures in the MS
fig3<-ggarrange(RP_bias, RP_CV, common.legend = TRUE, ncol=1, nrow=2, labels = c("a", "b"))
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig3.png"), fig3, width = 6, height = 5)
fig3
Indifference Pointtemp_data = sumexpSub_diff %>%
Indifference Point(linear regression)
#Observed Indifference Point
obs_model <- function(df) {
lm(repDur ~ curDur, data = df)
}
obs_InP <- AllExpData %>% #filter(curDur %in% c(1.1,1.4)) %>%
dplyr::group_by(NSub, Exp, WMSize) %>% nest() %>%
mutate(model = map(data, obs_model)) %>% # linear regression
mutate(slope = map(model, broom::tidy)) %>% # get estimates
unnest(slope, .drop = TRUE) %>% # remove raw data
select(-std.error,-statistic, -p.value) %>% # remove unnessary columns
spread(term, estimate) %>% # spread stimates
dplyr::rename(Intercept = `(Intercept)`, slope = curDur) # rename columns
obs_InP$model = NULL
obs_InP$data = NULL
obs_InP$inP = obs_InP$Intercept /(1-obs_InP$slope)
write.csv(obs_InP, file = paste0(modelPath, "/rlt/obs_InP.csv"))
mobs_InP = obs_InP %>% group_by(Exp, WMSize)%>%
dplyr::summarise(n=n(),
m_Intercept = mean(Intercept),
se_Intercept= sd(Intercept)/sqrt(n-1),
m_inP = mean(inP),
se_inP = sd(inP)/sqrt(n-1),
m_slope = mean(slope),
se_slope = sd(slope)/sqrt(n-1))
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
plt_InP_linear<- ggplot(data = mobs_InP, aes(x= Exp, y=m_inP, color = WMSize))+
geom_line(stat = "identity",position = position_dodge(width = 0.2))+
geom_point(stat = "identity",position = position_dodge(width = 0.2))+
geom_errorbar(width=.3, aes(ymin = m_inP - se_inP, ymax = m_inP + se_inP), position = position_dodge(width = 0.2)) +theme_new+
labs(colour = "Memory Load")+colorSet3+
xlab(' ')+ylab("indifference point (s)")+guides(shape="none")+
theme(legend.position = "top")
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_InP_linear.png"), plt_InP_linear, width = 5, height = 5)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
plt_InP_linear
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
RP_Exp2 <- ggplot(data = m_predY%>% filter(Exp == 'Exp2'), aes(x = curDur, y = m_m_repDur, color=WMSize, shape = as.factor('Observation'))) +
geom_point(size=2, alpha = 0.5)+
geom_segment(data =mobs_InP%>%filter(Exp == 'Exp2'), aes(x = 0.3, y = 0.3*m_slope+m_Intercept, xend = 2.2, yend = 2.2*m_slope+m_Intercept), arrow = NULL)+ ##dotted lines
geom_point(data =mobs_InP%>% filter(Exp == 'Exp2'), aes(x = m_inP, y = 0.1, color=WMSize), shape=21) + ## Intersection
geom_segment(data =mobs_InP%>% filter(Exp == 'Exp2'), aes(x=m_inP, y = 0.1, xend = m_inP, yend = m_inP), arrow = NULL)+ ##vertical lines for Intersection
#geom_line(data= m_newY%>% filter(Exp == 'Exp2'), aes(x=curDur, y=m_mu_r, color=WMSize)) +
geom_abline(slope=1, intercept=0, linetype = 2)+
labs(x="Sample intervals (s)", y="Reproduction (s)", shape=" ", color = "Memory Load")+
theme_new+colorSet3+guides(shape="none")+ ylim(0,2)
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_Exp2.png"), RP_Exp2, width = 4, height = 4)
RP_Exp2
Indifference Point (bootstraps)
#### Estimated Indifference Point
mRep_model <- function(df) {
lm(mu_r ~ curDur, data = df)
}
#Observed Indifference Point
obs_model <- function(df) {
lm(repDur ~ curDur, data = df)
}
# Custom function to find predicted indifference point
getPredInP_boot <- function(df, idx){
vars <- c('NSub', 'Exp', 'WMSize')
gp_vars = syms(vars)
slopes <- df[idx, ] %>%
dplyr::group_by(!!!gp_vars) %>% nest() %>% # nested data
mutate(model = map(data, mRep_model)) %>% # linear regression
mutate(slope = map(model, broom::tidy)) %>% # get estimates out
unnest(slope, .drop = TRUE) %>% # remove raw data
select(-std.error,-statistic, -p.value) %>% # remove unnessary clumns
spread(term, estimate) %>% # spread stimates
dplyr::rename(minRP = `(Intercept)`, slope = curDur) # rename columns
slopes$inP = slopes$minRP /(1-slopes$slope)
return(slopes$inP)
}
# Custom function to find observed indifference point
getRPInP_boot <- function(df, idx){
vars <- c('NSub', 'Exp', 'WMSize')
gp_vars = syms(vars)
slopes <- df[idx, ] %>%
dplyr::group_by(!!!gp_vars) %>% nest() %>% # nested data
mutate(model = map(data, obs_model)) %>% # linear regression
mutate(slope = map(model, broom::tidy)) %>% # get estimates out
unnest(slope, .drop = TRUE) %>% # remove raw data
select(-std.error,-statistic, -p.value) %>% # remove unnessary clumns
spread(term, estimate) %>% # spread stimates
dplyr::rename(minRP = `(Intercept)`, slope = curDur) # rename columns
slopes$inP = slopes$minRP /(1-slopes$slope)
return(slopes$inP)
}
#calculate the bootstrapped 95% confidence intervals
generateCI = FALSE # tag for generation CI
if(generateCI){
cilist <- NULL
for(expname in unique(AllDat_predY$Exp)){
for(nsub in unique(AllDat_predY$NSub)){
for(WMsize in unique(AllDat_predY$WMSize)){
dat = AllDat_predY %>% filter(Exp == expname, NSub == nsub, WMSize ==WMsize)
set.seed(100)
num = 1000
bs_predInP <- boot(dat, getPredInP_boot, R = num)
bs_RPInP <- boot(dat, getRPInP_boot, R = num)
ci = data.frame(
Exp = expname,
NSub = nsub,
WMSize = WMsize,
sd_predInP_boot = sd(bs_predInP$t),
m_predInP_boot = median(bs_predInP$t),
sd_RPInP_boot = sd(bs_RPInP$t),
mRP_InP_boot = median(bs_RPInP$t)
)
cilist = data.frame(rbind(cilist, ci))
}
}
}
write.csv(cilist, file = paste0("ci_list_median_1000.csv"))
}
# load the generated indifference point values and mark the outlier
cilist = read.csv(paste0("ci_list_median_1000.csv"))
cilist$Exp = as.factor(cilist$Exp)
cilist$WMSize= as.factor(cilist$WMSize)
cilist$inPOutlier = FALSE
cilist[which(cilist$mRP_InP_boot > 1.7 |cilist$mRP_InP_boot < 0.5 | cilist$m_predInP_boot < 0.5|cilist$m_predInP_boot > 1.7),"inPOutlier"] = TRUE
#check if the outlier is the same as the outliers in variable slope_pr
cilist%>% filter(inPOutlier == TRUE)
Inifference Point (curve)
AllDat_newY$predErr = AllDat_newY$mu_r - AllDat_newY$curDur
temp_newY <- AllDat_newY %>% filter(curDur > 0.8, curDur < 1.1) %>% select(Exp, WMSize, NSub, predErr, curDur)
InP_curve<- temp_newY%>% dplyr::group_by(Exp, WMSize, NSub)%>%
dplyr::summarise(minErr = min(abs(predErr)), idx = which.min(abs(predErr)))
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
InP_curve$InP_curve = temp_newY[InP_curve$idx,]$curDur
InP_curve$y = temp_newY[InP_curve$idx,]$predErr + temp_newY[InP_curve$idx,]$curDur
InP_curve
Inifference Point (Equation)
#sig_sm2 = sig_s2 + ls*(size[i]);
#wp_pr[i] = sig_sm2 /(sig_sm2 + sigma_pr2);
Bayparlist$wp_1 = (Bayparlist$sig_s2 + Bayparlist$ls *1 )/ (Bayparlist$sig_s2 + Bayparlist$ls *1 + Bayparlist$sig_pr2_log)
Bayparlist$wp_3 = (Bayparlist$sig_s2 + Bayparlist$ls *2) / (Bayparlist$sig_s2 + Bayparlist$ls *2 + Bayparlist$sig_pr2_log)
Bayparlist$wp_5 = (Bayparlist$sig_s2 + Bayparlist$ls *3) / (Bayparlist$sig_s2 + Bayparlist$ls *3 + Bayparlist$sig_pr2_log)
#log(IP)=[kr*M + (1-w_p)ks*M+ w_p*mu_p]/w_p
size = 1
Bayparlist$InP_1 = exp((Bayparlist$kr*size + (1-Bayparlist$wp_1)*Bayparlist$ks*size + Bayparlist$wp_1* Bayparlist$mu_pr_log)/Bayparlist$wp_1)
size = 2
Bayparlist$InP_3 = exp((Bayparlist$kr*size + (1-Bayparlist$wp_3)*Bayparlist$ks*size + Bayparlist$wp_3* Bayparlist$mu_pr_log)/Bayparlist$wp_3)
size = 3
Bayparlist$InP_5 = exp((Bayparlist$kr*size + (1-Bayparlist$wp_5)*Bayparlist$ks*size + Bayparlist$wp_5* Bayparlist$mu_pr_log)/Bayparlist$wp_5)
Bayparlist$sig2_post_1 = (Bayparlist$sig_pr2_log *(Bayparlist$sig_s2 +Bayparlist$ls*1))/(Bayparlist$sig_pr2_log + Bayparlist$sig_s2 +Bayparlist$ls*1) *0.5
Bayparlist$sig2_post_3 =(Bayparlist$sig_pr2_log *(Bayparlist$sig_s2 +Bayparlist$ls*3))/(Bayparlist$sig_pr2_log + Bayparlist$sig_s2 +Bayparlist$ls*3) *0.5
Bayparlist$sig2_post_5 = (Bayparlist$sig_pr2_log *(Bayparlist$sig_s2 +Bayparlist$ls*5))/(Bayparlist$sig_pr2_log + Bayparlist$sig_s2 +Bayparlist$ls*5) *0.5
Bayparlist%>%filter(Exp == 'Exp2')%>%select("NSub", "Exp", "sig2_post_1", "sig2_post_3","sig2_post_5")
Bayparlist$InP_old_1 = exp((Bayparlist$kr*size - (1-Bayparlist$wp_1)*Bayparlist$ks*size +(Bayparlist$sig_pr2_log *(Bayparlist$sig_s2 +Bayparlist$ls*1))/(Bayparlist$sig_pr2_log + Bayparlist$sig_s2 +Bayparlist$ls*1) *0.5 + Bayparlist$wp_1* Bayparlist$mu_pr_log)/Bayparlist$wp_1)
Bayparlist$InP_old_3 = exp((Bayparlist$kr*size - (1-Bayparlist$wp_3)*Bayparlist$ks*size +(Bayparlist$sig_pr2_log *(Bayparlist$sig_s2 +Bayparlist$ls*3))/(Bayparlist$sig_pr2_log + Bayparlist$sig_s2 +Bayparlist$ls*3) *0.5 + Bayparlist$wp_3* Bayparlist$mu_pr_log)/Bayparlist$wp_3)
Bayparlist$InP_old_5 = exp((Bayparlist$kr*size - (1-Bayparlist$wp_5)*Bayparlist$ks*size +(Bayparlist$sig_pr2_log *(Bayparlist$sig_s2 +Bayparlist$ls*5))/(Bayparlist$sig_pr2_log + Bayparlist$sig_s2 +Bayparlist$ls*5) *0.5 + Bayparlist$wp_5* Bayparlist$mu_pr_log)/Bayparlist$wp_5)
Bayparlist%>%select(NSub, Exp, wp_1, wp_3, wp_5, InP_1, InP_3, InP_5, InP_old_1, InP_old_3, InP_old_5)%>%
unite(newcol1, wp_1, InP_1, InP_old_1)%>%
unite(newcol3, wp_3, InP_3, InP_old_3)%>%
unite(newcol5, wp_5, InP_5, InP_old_5)%>%
melt(id.vars = c("NSub", "Exp")) %>%
dplyr::rename(
WMSize = variable,
newcol = value
)%>%
separate(newcol, c("wp", "InP_Eq", "InP_Eq_old"), sep = "_") -> Baypar_WMSize
Baypar_WMSize$WMSize <- substr(Baypar_WMSize$WMSize, 7, 7)
Baypar_WMSize$WMSize <- as.factor(Baypar_WMSize$WMSize)
Baypar_WMSize$InP_Eq = as.numeric(Baypar_WMSize$InP_Eq)
Baypar_WMSize$InP_Eq_old = as.numeric(Baypar_WMSize$InP_Eq_old)
levels(Baypar_WMSize$WMSize) = c("low", "medium", "high")
##save csv for SPSS
left_join(Baypar_WMSize, InP_curve%>%select("Exp", "WMSize","NSub", "InP_curve"), by = c("Exp","WMSize","NSub")) -> Baypar_WMSize
left_join(Baypar_WMSize, cilist %>% select("Exp", "NSub", "WMSize", "m_predInP_boot", "mRP_InP_boot"), by = c("Exp","WMSize","NSub"))-> Baypar_WMSize
check indifference Points
Baypar_WMSize_melt <- melt(Baypar_WMSize%>%select("Exp","WMSize","NSub", "InP_Eq", "InP_Eq_old", "InP_curve", "m_predInP_boot", "mRP_InP_boot"), id.vars = c("Exp","WMSize","NSub"),
variable.name = "Type",
value.name = "InP" )
Baypar_WMSize_melt %>%filter(Type== "InP_curve") %>% dplyr::group_by(Exp)%>%
dplyr::summarise(n = n(),
m_InP = mean(InP), se_InP= sd(InP)/sqrt(n-1))
Baypar_WMSize_melt$Exp = as.factor(Baypar_WMSize_melt$Exp)
Baypar_WMSize_melt%>% dplyr::group_by(Exp, WMSize, Type)%>%
dplyr::summarise(n = n(),
m_InP = mean(InP), se_InP= sd(InP)/sqrt(n-1))-> mBaypar_InP
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
#plot InP_curve(the intersections of the Prediction curve with the diagonal)
plt_InP<- ggplot(data = mBaypar_InP%>%filter(Type== "InP_curve"), aes(x= Exp, y=m_InP, color = WMSize, shape = Type))+
geom_line(stat = "identity",position = position_dodge(width = 0.2))+
geom_point(stat = "identity",position = position_dodge(width = 0.2))+
geom_errorbar(width=.2, aes(ymin = m_InP - se_InP, ymax = m_InP + se_InP), position = position_dodge(width = 0.2)) +theme_new+
labs(colour = "Memory Load")+colorSet3+
xlab(' ')+ylab("indifference point (s)")+guides(shape="none")+
theme(legend.position = "top")
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_InP.png"), plt_InP, width = 3, height = 3)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
plt_InP
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
temp_inp = mBaypar_InP%>%filter(Type== "InP_curve", Exp == 'Exp2')
RP_bias_Exp2 <- ggplot(data = m_predY%>%filter(Exp == 'Exp2'), aes(x = curDur, y = m_m_repDur - curDur,color=WMSize, shape = as.factor('Observation'))) +
geom_point(size=2, alpha = 0.5)+
geom_line(data= m_newY%>%filter(Exp == 'Exp2'), aes(x=curDur, y=m_mu_r-curDur, color=WMSize)) +
geom_hline(yintercept = 0, linetype='dashed')+
geom_point(data =temp_inp, aes(x = m_InP+0.03, y = -0.5, color=WMSize), shape=21) + ## Intersection
geom_segment(data =temp_inp, aes(x=m_InP+0.03, y = -0.5, xend = m_InP+0.03, yend = 0), arrow = NULL)+ ##vertical lines for Intersection
labs(x=" ", y="Reproduction bias (s)", shape=" ", color = "Memory Load")+theme_new+
colorSet3+guides(shape="none")
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_bias_Exp2.png"), RP_bias_Exp2, width = 5, height = 3)
RP_bias_Exp2
anova on observed InP
ezANOVA(data = Baypar_WMSize, dv= InP_curve, wid=NSub, within = .(WMSize), between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 Exp 4 75 1.198179 3.187005e-01 0.054696877
## 3 WMSize 2 150 2.377013 9.631587e-02 0.002987228
## 4 Exp:WMSize 8 150 4.618376 4.655267e-05 * 0.022755616
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 3 WMSize 0.3471828 1.001507e-17 *
## 4 Exp:WMSize 0.3471828 1.001507e-17 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 3 WMSize 0.6050276 0.1210006309 0.6095966 0.1207130064
## 4 Exp:WMSize 0.6050276 0.0009604633 * 0.6095966 0.0009269622 *
pairwise.t.test(Baypar_WMSize$InP_curve, Baypar_WMSize$Exp)
##
## Pairwise comparisons using t tests with pooled SD
##
## data: Baypar_WMSize$InP_curve and Baypar_WMSize$Exp
##
## Exp1 Exp2 Exp3 Exp4a
## Exp2 1.000 - - -
## Exp3 1.000 1.000 - -
## Exp4a 1.000 0.951 1.000 -
## Exp4b 0.075 0.019 0.031 0.604
##
## P value adjustment method: holm
ezANOVA(data = Baypar_WMSize %>% filter(Exp =='Exp2'), dv= InP_curve, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 WMSize 2 30 5.834043 0.007240452 * 0.0340583
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 WMSize 0.1674321 3.688675e-06 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 WMSize 0.5456824 0.02548011 * 0.5558358 0.02476867 *
ezANOVA(data = Baypar_WMSize %>% filter(Exp =='Exp3'), dv= InP_curve, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 WMSize 2 30 8.753035 0.001012937 * 0.03895252
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 WMSize 0.06186308 3.467534e-09 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 WMSize 0.5159594 0.009073344 * 0.5194236 0.008930414 *
ezANOVA(data = Baypar_WMSize %>% filter(Exp =='Exp4a'), dv= InP_curve, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 WMSize 2 30 0.7386935 0.4862266 0.00420096
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 WMSize 0.07926534 1.965999e-08 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 WMSize 0.5206341 0.408452 0.5251298 0.4094818
ezANOVA(data = Baypar_WMSize %>% filter(Exp =='Exp4b'), dv= InP_curve, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 WMSize 2 30 3.57389 0.04053459 * 0.04484435
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 WMSize 0.450539 0.003768136 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 WMSize 0.6453857 0.0647428 0.6802842 0.06184001
Export data for spss
### change the table struction for spss
Bayparlist_list <-split(Baypar_WMSize %>% select(c("NSub","Exp","WMSize", "wp","InP_curve")), Baypar_WMSize$WMSize)
Bayparlist_spss = NULL
for (i in 1: length(Bayparlist_list)){
temp = Bayparlist_list[[i]]
WMSize = unique(temp$WMSize)
temp$WMSize = NULL
colnames(temp) = c('NSub','Exp', paste0('wp_',WMSize), paste0('InP_curve_',WMSize))
if(i == 1)
Bayparlist_spss = temp
else
Bayparlist_spss = left_join(Bayparlist_spss, temp, by=c("Exp", "NSub"))
}
write_csv(Bayparlist_spss, paste0(modelPath, '/rlt/Bayparlist_spss.csv'))
anova analysis
Anova on mean reproduction biases
mpredY_sub$RP_Bias = mpredY_sub$m_repDur-mpredY_sub$curDur
RP_bias_Anova <- ezANOVA(data = mpredY_sub, dv= RP_Bias, wid=NSub, within= .(curDur, WMSize), between = .(Exp), detailed = TRUE, return_aov = FALSE )
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: "curDur" will be treated as numeric.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Warning: There is at least one numeric within variable, therefore aov() will be
## used for computation and no assumption checks will be obtained.
RP_bias_Anova
## $ANOVA
## Effect DFn DFd SSn SSd F p p<.05
## 1 Exp 4 75 0.37301817 6.5861689 1.061936 3.813811e-01
## 2 curDur 1 75 49.69246885 7.6805406 485.243856 1.753714e-34 *
## 4 WMSize 2 150 0.04636729 0.6908236 5.033914 7.656681e-03 *
## 3 Exp:curDur 4 75 0.78848535 7.6805406 1.924878 1.150304e-01
## 5 Exp:WMSize 8 150 0.25333639 0.6908236 6.875933 1.093463e-07 *
## 6 curDur:WMSize 2 150 0.12351691 0.3507396 26.412100 1.488886e-10 *
## 7 Exp:curDur:WMSize 8 150 0.12715860 0.3507396 6.797704 1.342826e-07 *
## ges
## 1 0.023787466
## 2 0.764490798
## 4 0.003019758
## 3 0.048984109
## 5 0.016279575
## 6 0.008004056
## 7 0.008238098
# main effect of Duration F(1.177, 3.532) = 377.965, p < .001, ηp² = .863.
(RP_bias_Anova$ANOVA)$DFn[3] *(RP_bias_Anova$`Sphericity Corrections`)$GGe[1]
## numeric(0)
(RP_bias_Anova$ANOVA)$DFd[3] *(RP_bias_Anova$`Sphericity Corrections`)$GGe[1]
## numeric(0)
#Duration × Experiment, F(12, 240) = 2.506, p = .004, ηp² = .111
(RP_bias_Anova$ANOVA)$DFn[5] *(RP_bias_Anova$`Sphericity Corrections`)$GGe[2]
## numeric(0)
(RP_bias_Anova$ANOVA)$DFd[5] *(RP_bias_Anova$`Sphericity Corrections`)$GGe[2]
## numeric(0)
mpredY_sub <- ungroup(mpredY_sub)
res.aov <- rstatix::anova_test(data = mpredY_sub, dv = RP_Bias, wid = NSub, within = c(curDur, WMSize), between = Exp)
## Warning: The 'wid' column contains duplicate ids across between-subjects
## variables. Automatic unique id will be created
get_anova_table(res.aov, correction = "GG")
variance of prior (anova)
ezANOVA(data = Bayparlist, dv= sig_pr2_log, wid=NSub, between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Coefficient covariances computed by hccm()
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 1 Exp 4 75 0.7160159 0.5836025 0.03678287
##
## $`Levene's Test for Homogeneity of Variance`
## DFn DFd SSn SSd F p p<.05
## 1 4 75 0.1582023 3.779256 0.7848881 0.5385686
variance of motor noise (anova)
ezANOVA(data = Bayparlist, dv= sig_mn2, wid=NSub, between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Coefficient covariances computed by hccm()
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 1 Exp 4 75 1.864222 0.1255674 0.09043378
##
## $`Levene's Test for Homogeneity of Variance`
## DFn DFd SSn SSd F p p<.05
## 1 4 75 0.0005070715 0.007658422 1.241456 0.3007183
Ks anova
ezANOVA(data = Bayparlist, dv= ks, wid=NSub, between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Coefficient covariances computed by hccm()
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 1 Exp 4 75 21.71153 6.283521e-12 * 0.5365969
##
## $`Levene's Test for Homogeneity of Variance`
## DFn DFd SSn SSd F p p<.05
## 1 4 75 0.4637314 1.347294 6.453651 0.0001614605 *
mean of prior (anova)
ezANOVA(data = Bayparlist, dv= mu_pr, wid=NSub, between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Coefficient covariances computed by hccm()
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 1 Exp 4 75 1.372231 0.2516599 0.06819476
##
## $`Levene's Test for Homogeneity of Variance`
## DFn DFd SSn SSd F p p<.05
## 1 4 75 0.8682032 27.17694 0.5989934 0.6644974
weight of prior
Weight of the prior \(w_p\)
Baypar_WMSize$wp = as.numeric(Baypar_WMSize$wp)
mwp <- Baypar_WMSize%>%dplyr::group_by(Exp, WMSize)%>% dplyr::summarise(m_wp = mean(wp), n= n(), se_wp = sd(wp)/sqrt(n-1))
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
plt_wp <- ggplot(mwp, aes(Exp, m_wp, ymin = m_wp - se_wp, ymax = m_wp + se_wp, group =interaction(Exp, WMSize), color = factor(WMSize)), shape = factor(WMSize))+
geom_line(stat = "identity",position = position_dodge(width = 0.2))+
geom_point(stat = "identity",position = position_dodge(width = 0.2))+
geom_errorbar(width=.2, position = position_dodge(width = 0.2)) +
#coord_cartesian(ylim = c(0.5, 1)) +
colorSet5+
labs(x = "", y = TeX("Weight of the prior $w_p$"), color = 'Memory Load') +
theme_new + theme(legend.position="top")
plt_wp
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_wp.png"), plt_wp, width = 3, height = 3)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
combine InP and wp
fig4<-ggarrange(plt_wp, plt_InP, common.legend = TRUE, ncol=2, nrow=1, labels = c("a", "b"))
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig4.png"), fig4, width = 6, height = 3)
fig4
fig5<-ggarrange(RP_bias_obs, plt_InP_linear, common.legend = TRUE, ncol=2, nrow=1, widths = c(5,3), labels = c("a", "b"))
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig5.png"), fig5, width = 6, height = 3)
fig5
ANOVA analysis on wp
ezANOVA(data = Baypar_WMSize, dv= wp, wid=NSub, within=.(WMSize), between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 Exp 4 75 2.249548 7.163591e-02 0.106466774
## 3 WMSize 2 150 123.376009 2.079259e-32 * 0.011162610
## 4 Exp:WMSize 8 150 25.074955 3.255859e-24 * 0.009093747
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 3 WMSize 0.0228096 1.779717e-61 *
## 4 Exp:WMSize 0.0228096 1.779717e-61 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 3 WMSize 0.5057682 1.120300e-17 * 0.506003 1.102351e-17 *
## 4 Exp:WMSize 0.5057682 2.492781e-13 * 0.506003 2.463140e-13 *
ezANOVA(data = Baypar_WMSize%>%filter(Exp =='Exp2'), dv= wp, wid=NSub, within= .(WMSize) )
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 WMSize 2 30 55.0847 9.057763e-11 * 0.04339272
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 WMSize 0.02641187 8.965915e-12 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 WMSize 0.5066913 1.869688e-06 * 0.508133 1.815792e-06 *
ezANOVA(data = Baypar_WMSize%>%filter(Exp =='Exp4a'), dv= wp, wid=NSub, within= .(WMSize) )
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 WMSize 2 30 33.98248 1.953251e-08 * 0.01772593
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 WMSize 0.01054139 1.44639e-14 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 WMSize 0.5026493 3.187613e-05 * 0.5032182 3.160504e-05 *
ezANOVA(data = Baypar_WMSize%>%filter(Exp =='Exp4b'), dv= wp, wid=NSub, within= .(WMSize) )
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 WMSize 2 30 34.92081 1.469415e-08 * 0.02125106
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 WMSize 0.02260196 3.013161e-12 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 WMSize 0.5057151 2.627002e-05 * 0.5069454 2.578052e-05 *
### pairwise.t.test on wp
pairwise.t.test(Baypar_WMSize$wp, Baypar_WMSize$Exp)
##
## Pairwise comparisons using t tests with pooled SD
##
## data: Baypar_WMSize$wp and Baypar_WMSize$Exp
##
## Exp1 Exp2 Exp3 Exp4a
## Exp2 0.00016 - - -
## Exp3 0.12532 0.21889 - -
## Exp4a 1.00000 0.00029 0.16667 -
## Exp4b 1.00000 0.00323 0.47867 1.00000
##
## P value adjustment method: holm
##independent T test
t.test((Baypar_WMSize%>%filter(Exp %in%c('Exp1')))$wp)
##
## One Sample t-test
##
## data: (Baypar_WMSize %>% filter(Exp %in% c("Exp1")))$wp
## t = 14.032, df = 47, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.3459136 0.4617008
## sample estimates:
## mean of x
## 0.4038072
t.test((Baypar_WMSize%>%filter(Exp %in%c('Exp2')))$wp)
##
## One Sample t-test
##
## data: (Baypar_WMSize %>% filter(Exp %in% c("Exp2")))$wp
## t = 16.632, df = 47, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.5184761 0.6611631
## sample estimates:
## mean of x
## 0.5898196
t.test((Baypar_WMSize%>%filter(Exp %in%c('Exp3')))$wp)
##
## One Sample t-test
##
## data: (Baypar_WMSize %>% filter(Exp %in% c("Exp3")))$wp
## t = 20.765, df = 47, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.4554826 0.5532064
## sample estimates:
## mean of x
## 0.5043445
t.test((Baypar_WMSize%>%filter(Exp %in%c('Exp4a')))$wp)
##
## One Sample t-test
##
## data: (Baypar_WMSize %>% filter(Exp %in% c("Exp4a")))$wp
## t = 14.029, df = 47, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.3520475 0.4699193
## sample estimates:
## mean of x
## 0.4109834
t.test((Baypar_WMSize%>%filter(Exp %in%c('Exp4b')))$wp)
##
## One Sample t-test
##
## data: (Baypar_WMSize %>% filter(Exp %in% c("Exp4b")))$wp
## t = 14.529, df = 47, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.3777757 0.4992013
## sample estimates:
## mean of x
## 0.4384885
standard variance of Ds \(\sigma_{s}\)
ezANOVA(data = Bayparlist, dv= sig_s2, wid=NSub, between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Coefficient covariances computed by hccm()
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 1 Exp 4 75 0.3845193 0.8190605 0.02009558
##
## $`Levene's Test for Homogeneity of Variance`
## DFn DFd SSn SSd F p p<.05
## 1 4 75 0.02897922 0.8327796 0.6524659 0.6269356
### pairwise.t.test on standard variance of $D_s$
pairwise.t.test(Bayparlist$sig_s2, Bayparlist$Exp)
##
## Pairwise comparisons using t tests with pooled SD
##
## data: Bayparlist$sig_s2 and Bayparlist$Exp
##
## Exp1 Exp2 Exp3 Exp4a
## Exp2 1 - - -
## Exp3 1 1 - -
## Exp4a 1 1 1 -
## Exp4b 1 1 1 1
##
## P value adjustment method: holm
##independent T test
t.test((Bayparlist%>%filter(Exp %in%c('Exp1')))$sig_s2)
##
## One Sample t-test
##
## data: (Bayparlist %>% filter(Exp %in% c("Exp1")))$sig_s2
## t = 7.2797, df = 15, p-value = 2.698e-06
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.04681825 0.08558531
## sample estimates:
## mean of x
## 0.06620178
t.test((Bayparlist%>%filter(Exp %in%c('Exp2')))$sig_s2)
##
## One Sample t-test
##
## data: (Bayparlist %>% filter(Exp %in% c("Exp2")))$sig_s2
## t = 2.5569, df = 15, p-value = 0.0219
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.01307008 0.14402351
## sample estimates:
## mean of x
## 0.07854679
t.test((Bayparlist%>%filter(Exp %in%c('Exp3')))$sig_s2)
##
## One Sample t-test
##
## data: (Bayparlist %>% filter(Exp %in% c("Exp3")))$sig_s2
## t = 11.647, df = 15, p-value = 6.497e-09
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.06165666 0.08927951
## sample estimates:
## mean of x
## 0.07546809
t.test((Bayparlist%>%filter(Exp %in%c('Exp4a')))$sig_s2)
##
## One Sample t-test
##
## data: (Bayparlist %>% filter(Exp %in% c("Exp4a")))$sig_s2
## t = 7.0467, df = 15, p-value = 3.959e-06
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.03222067 0.06016474
## sample estimates:
## mean of x
## 0.04619271
t.test((Bayparlist%>%filter(Exp %in%c('Exp4b')))$sig_s2)
##
## One Sample t-test
##
## data: (Bayparlist %>% filter(Exp %in% c("Exp4b")))$sig_s2
## t = 1.7736, df = 15, p-value = 0.09643
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.01870072 0.20404884
## sample estimates:
## mean of x
## 0.09267406
# calculate conhens d
cohensD(sig_s2 ~ Exp,
data = Bayparlist%>%filter(Exp %in% c('Exp1', 'Exp2')),
method = "paired")
## Warning in cohensD(sig_s2 ~ Exp, data = Bayparlist %>% filter(Exp %in%
## c("Exp1", : calculating paired samples Cohen's d using formula input. Results
## will be incorrect if cases do not appear in the same order for both levels of
## the grouping factor
## [1] 0.08918868
cohensD(sig_s2 ~ Exp,
data = Bayparlist%>%filter(Exp %in% c('Exp2', 'Exp3')),
method = "paired")
## Warning in cohensD(sig_s2 ~ Exp, data = Bayparlist %>% filter(Exp %in%
## c("Exp2", : calculating paired samples Cohen's d using formula input. Results
## will be incorrect if cases do not appear in the same order for both levels of
## the grouping factor
## [1] 0.02474716
cohensD(sig_s2 ~ Exp,
data = Bayparlist%>%filter(Exp %in% c('Exp3', 'Exp4a')),
method = "paired")
## Warning in cohensD(sig_s2 ~ Exp, data = Bayparlist %>% filter(Exp %in%
## c("Exp3", : calculating paired samples Cohen's d using formula input. Results
## will be incorrect if cases do not appear in the same order for both levels of
## the grouping factor
## [1] 0.8774215
cohensD(sig_s2 ~ Exp,
data = Bayparlist%>%filter(Exp %in% c('Exp2', 'Exp4a')),
method = "paired")
## Warning in cohensD(sig_s2 ~ Exp, data = Bayparlist %>% filter(Exp %in%
## c("Exp2", : calculating paired samples Cohen's d using formula input. Results
## will be incorrect if cases do not appear in the same order for both levels of
## the grouping factor
## [1] 0.2541347
cv slope
ezANOVA(data = CVslopes, dv= obs_cv_slope, wid=NSub, between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Warning: Collapsing data to cell means. *IF* the requested effects are a subset
## of the full design, you must use the "within_full" argument, else results may be
## inaccurate.
## Coefficient covariances computed by hccm()
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 1 Exp 4 75 0.08928208 0.9855437 0.004739144
##
## $`Levene's Test for Homogeneity of Variance`
## DFn DFd SSn SSd F p p<.05
## 1 4 75 0.08788252 1.346571 1.223699 0.3079875
### pairwise.t.test on standard variance of $D_s$
pairwise.t.test(CVslopes$obs_cv_slope, CVslopes$Exp)
##
## Pairwise comparisons using t tests with pooled SD
##
## data: CVslopes$obs_cv_slope and CVslopes$Exp
##
## Exp1 Exp2 Exp3 Exp4a
## Exp2 1 - - -
## Exp3 1 1 - -
## Exp4a 1 1 1 -
## Exp4b 1 1 1 1
##
## P value adjustment method: holm
##independent T test
t.test((CVslopes%>%filter(Exp %in%c('Exp1')))$obs_cv_slope)
##
## One Sample t-test
##
## data: (CVslopes %>% filter(Exp %in% c("Exp1")))$obs_cv_slope
## t = -6.5917, df = 47, p-value = 3.406e-08
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.2970793 -0.1581479
## sample estimates:
## mean of x
## -0.2276136
t.test((CVslopes%>%filter(Exp %in%c('Exp2')))$obs_cv_slope)
##
## One Sample t-test
##
## data: (CVslopes %>% filter(Exp %in% c("Exp2")))$obs_cv_slope
## t = -7.47, df = 47, p-value = 1.59e-09
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.3369816 -0.1939864
## sample estimates:
## mean of x
## -0.265484
t.test((CVslopes%>%filter(Exp %in%c('Exp3')))$obs_cv_slope)
##
## One Sample t-test
##
## data: (CVslopes %>% filter(Exp %in% c("Exp3")))$obs_cv_slope
## t = -6.932, df = 47, p-value = 1.037e-08
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.3153826 -0.1735021
## sample estimates:
## mean of x
## -0.2444424
t.test((CVslopes%>%filter(Exp %in%c('Exp4a')))$obs_cv_slope)
##
## One Sample t-test
##
## data: (CVslopes %>% filter(Exp %in% c("Exp4a")))$obs_cv_slope
## t = -5.4143, df = 47, p-value = 2.047e-06
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.3547829 -0.1625578
## sample estimates:
## mean of x
## -0.2586703
t.test((CVslopes%>%filter(Exp %in%c('Exp4b')))$obs_cv_slope)
##
## One Sample t-test
##
## data: (CVslopes %>% filter(Exp %in% c("Exp4b")))$obs_cv_slope
## t = -7.1687, df = 47, p-value = 4.537e-09
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.2997242 -0.1683650
## sample estimates:
## mean of x
## -0.2340446
# calculate conhens d
cohensD(obs_cv_slope ~ Exp,
data = CVslopes%>%filter(Exp %in% c('Exp1', 'Exp2')),
method = "paired")
## Warning in cohensD(obs_cv_slope ~ Exp, data = CVslopes %>% filter(Exp %in% :
## calculating paired samples Cohen's d using formula input. Results will be
## incorrect if cases do not appear in the same order for both levels of the
## grouping factor
## [1] 0.1178123
cohensD(obs_cv_slope ~ Exp,
data = CVslopes%>%filter(Exp %in% c('Exp2', 'Exp3')),
method = "paired")
## Warning in cohensD(obs_cv_slope ~ Exp, data = CVslopes %>% filter(Exp %in% :
## calculating paired samples Cohen's d using formula input. Results will be
## incorrect if cases do not appear in the same order for both levels of the
## grouping factor
## [1] 0.06429713
cohensD(obs_cv_slope ~ Exp,
data = CVslopes%>%filter(Exp %in% c('Exp3', 'Exp4a')),
method = "paired")
## Warning in cohensD(obs_cv_slope ~ Exp, data = CVslopes %>% filter(Exp %in% :
## calculating paired samples Cohen's d using formula input. Results will be
## incorrect if cases do not appear in the same order for both levels of the
## grouping factor
## [1] 0.04045837
cohensD(obs_cv_slope ~ Exp,
data = CVslopes%>%filter(Exp %in% c('Exp2', 'Exp4a')),
method = "paired")
## Warning in cohensD(obs_cv_slope ~ Exp, data = CVslopes %>% filter(Exp %in% :
## calculating paired samples Cohen's d using formula input. Results will be
## incorrect if cases do not appear in the same order for both levels of the
## grouping factor
## [1] 0.01404826
Model prediction error
m_predErr_sub<- mpredY_sub%>%
dplyr::group_by(Exp, WMSize, NSub) %>% dplyr::summarise(
mpredRP_err=mean(predRP_err),
mpredVar_err=mean(predVar_err),
mpredcv_err = mean(predcv_err),
mpredRP_rerr = mean(predRP_rerr),
mpredVar_rerr = mean(predVar_rerr),
mpredcv_rerr = mean(predcv_rerr))
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
m_predErr<- m_predY%>%
dplyr::group_by(Exp, WMSize) %>% dplyr::summarise(
mmpredcv_err = mean(mpredcv_err),
mmpredRP_err=mean(mpredRP_err),
mmpredVar_err=mean(mpredVar_err),
mmpredRP_rerr = mean(mpredRP_rerr),
mmpredVar_rerr = mean(mpredVar_rerr),
mmpredcv_rerr = mean(mpredcv_rerr))
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
m_predErr
plt_rErrorScatter = ggplot(m_predErr_sub, aes(mpredRP_rerr*100, mpredcv_rerr*100, color = WMSize, alpha = .9)) +
#geom_hline(yintercept = 0, linetype='dashed')+ geom_vline(xintercept = 0, linetype='dashed')+
geom_point() +
geom_point(data = m_predErr, aes(mmpredRP_rerr*100, mmpredcv_rerr*100, color = WMSize, alpha = .9, size = 1 ))+
xlab('Relative prediction error of reproduction mean(%)')+ ylab('Relative prediction error of CV (%)')+colorSet3+
facet_wrap(~Exp)+
theme_new+ theme(legend.position = 'top')+guides(size="none")+guides(alpha="none")
plt_rErrorScatter
plt_rErrorScatter = ggplot(m_predErr_sub, aes(mpredRP_rerr*100, mpredVar_rerr*100, color = WMSize, alpha = .9)) +
#geom_hline(yintercept = 0, linetype='dashed')+ geom_vline(xintercept = 0, linetype='dashed')+
geom_point() +
geom_point(data = m_predErr, aes(mmpredRP_rerr*100, mmpredVar_rerr*100, color = WMSize, alpha = .9, size = 1 ))+
xlab('Relative prediction error in the RP means (%)')+ ylab('Relative prediction error in the RP variance (%)')+colorSet3+
facet_wrap(~Exp)+
theme_new+ theme(legend.position = 'top')+guides(size="none")+guides(alpha="none")
plt_rErrorScatter
plt_ErrorScatter = ggplot(m_predErr_sub, aes(mpredRP_err, mpredVar_err, color = WMSize, alpha = .9)) +
geom_hline(yintercept = 0, linetype='dashed')+ geom_vline(xintercept = 0, linetype='dashed')+
geom_point() +
geom_point(data = m_predErr, aes(mmpredRP_err, mmpredVar_err, color = WMSize, alpha = .9, size = 1 ))+
xlab('Prediction error in the RP means (s)')+ ylab('Prediction error in the RP variance (s)')+colorSet3+
facet_wrap(~Exp)+
theme_new+ theme(legend.position = 'top')+guides(size="none")+guides(alpha="none")
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_ErrorScatter.png"), plt_ErrorScatter, width = 7, height = 4)
plt_ErrorScatter
model comparison (logarithmic vs. linear)
m_predErr_sub$model = 'logarithmic'
m_predErr$model = 'logarithmic'
linear_model = 'linear_rstan'
m_predErr_linear = read.csv(paste0(getwd(), "/", rstanmodelPath, '/models/', linear_model, "/rlt/m_predErr_", linear_model, ".csv"))
m_predErr_linear$X = NULL
m_predErr_sub_linear = read.csv(paste0(getwd(), "/", rstanmodelPath, '/models/', linear_model, "/rlt/m_predErr_sub_", linear_model, ".csv"))
m_predErr_sub_linear$X = NULL
m_predErr_sub_all = rbind(m_predErr_sub, m_predErr_sub_linear)
m_predErr_all = rbind(m_predErr, m_predErr_linear)
m_predErr_all$WMSize = as.factor(m_predErr_all$WMSize)
levels(m_predErr_all$WMSize) = c("low", "medium", "high")
temp = m_predErr_all %>% filter(model == 'logarithmic') %>%summarise(abs_mmpredcv_err = abs(mmpredcv_err))
plt_Err_CV_all = ggplot(m_predErr_all, aes(abs(mmpredRP_err), abs(mmpredcv_err), group = interaction(model, Exp, WMSize), color = model, shape = Exp, size = WMSize)) +
geom_point() +
geom_hline(yintercept = round(max(temp$abs_mmpredcv_err), 4)+0.0005, linetype='dashed')+
xlab('Prediction error in the RP means (s)')+ ylab('Prediction error in CV')+colorSet3+
scale_shape_manual(values = c(6, 7, 16, 17,8)) +
theme_new+
theme(legend.position = 'top')+
labs(size = 'Memory Load')+
guides(colour = guide_legend(order = 1, nrow=2,byrow=TRUE),
shape = guide_legend(order =2, nrow=2,byrow=TRUE),
size = guide_legend(order = 3, nrow=3,byrow=TRUE))
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_Err_CV_all.png"), plt_Err_CV_all, width = 7, height = 4)
## Warning: Using size for a discrete variable is not advised.
plt_Err_CV_all
## Warning: Using size for a discrete variable is not advised.
m_predY_acc = m_predErr_sub_all%>%
dplyr::group_by(Exp, model) %>%
dplyr::summarize(mmpredRP_rerr = mean(mpredRP_rerr)*100,
mmpredVar_rerr = mean(mpredVar_rerr)*100,
mmpredcv_rerr = mean(mpredcv_rerr)*100,
mmpredRP_acc = (1-mean(mpredRP_rerr))*100,
mmpredVar_acc = (1-mean(mpredVar_rerr))*100,
mmpredCV_acc = (1-mean(mpredcv_rerr))*100)
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
m_predY_acc